function phaseout = GPOPSContinuous(input)

    settings = input.auxdata.settings;
%% counter for debugging
    if settings.counter
        load('counter.mat')
        counter = counter+1;
        save('counter.mat','counter');    
    end

    gpops_fun = input.auxdata.gpops_fun;
    var = gpops_fun.compute_var_from_input(input);

    %% Specify laws of motion
    phaseout.dynamics = gpops_fun.dynamics(var);
    
    %% integrals
    
    % welfare
    welfare_M = gpops_fun.welfare_participants_M(var) ./ input.auxdata.settings.inpscale.welfare;
    welfare_R = gpops_fun.welfare_participants_R(var) ./ input.auxdata.settings.inpscale.welfare;
    welfare_C = gpops_fun.welfare_participants_C(var) ./ input.auxdata.settings.inpscale.welfare;
    welfare_nonpart = gpops_fun.welfare_non_participants(var) ./ ...
        input.auxdata.settings.inpscale.welfare;
    welfare = gpops_fun.welfare(var) ./ input.auxdata.settings.inpscale.welfare;
    
    % Resource constraint
    resource_constraint = gpops_fun.resource_constraint(var);
    
    % Conistency conditions
    cc_M =  gpops_fun.labor_market_clearing_M(var);
    cc_R =  gpops_fun.labor_market_clearing_R(var);
    cc_C =  gpops_fun.labor_market_clearing_C(var);

    cons_M = gpops_fun.consumption_M(var);
    cons_R = gpops_fun.consumption_R(var);
    cons_C = gpops_fun.consumption_C(var);
    
    phaseout.integrand = [welfare_M, welfare_R, welfare_C, welfare_nonpart, resource_constraint, ...
                        cc_M, cc_R, cc_C, cons_M, cons_R, cons_C]; 

    
end
